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We investigate macroscopic quantum tunneling of a Bose condensate and how it is affected by the 
interatomic interaction. We study the dynamics of a condensate falling under gravity and scattering 
on a Gaussian potential barrier that models a mirror formed by a far-detuned sheet of light. We 
observe bouncing, interference and quantum tunneling of the condensate. We find that the tunneling 
rate is very sensitive to the interatomic interaction and to the shape of the condensate. Under many 
conditions the tunneling rate is strongly enhanced by the interaction as achieved, for instance, by 
increasing the number of condensed particles. In a quasi ID situation the tunneling pulse displays 
two peaks. The quantum tunneling can be quasi-periodic and in this way one could generate coherent 
Bose condensed atomic pulses. PACS numbers: 03.75.Fi, 32.80.Pj, 42.50.Vk 



I. INTRODUCTION 

Bose condensates, as achieved in trapped clouds of 
alkali- metal atoms [1-3], are ideal systems to study quan- 
tum phenomena at macroscopic level [4] , like interference 
[5] and diffraction of matter waves. Tunneling is one 
of the striking quantum phenomena, which is important 
also in many applications. In the case of non-interacting 
particles at very low temperatures, when almost all par- 
ticles are in the condensate, the tunneling of the con- 
densate has the same character as in the case of a sin- 
gle particle. In this paper we address the question of 
how tunneling is affected by the interatomic interaction. 
The case of a condensate in a double-well potential has 
been already theoretically investigated and it is known 
that the interatomic interaction tends to suppress tun- 
neling [6-10]. This effect can be so strong that under 
certain conditions it has been predicted that tunneling is 
suppressed completely and the condensate remains self- 
trapped in one of the minima of the potential. Here we 
show that under other conditions the tunneling rate can 
be enhanced by the interatomic interaction. 

The case we study is that of a condensate impinging 
on a potential barrier that could model a mirror formed 
by a far-detuned sheet of light. We find that the in- 
teratomic interaction and the geometrical aspects of the 
system, like the aspect ratio of the cloud or the fact the 
the cloud remains trapped in the transverse directions, 
have a strong effect on the tunneling probability. This 
can be either enhanced or depressed and this effect can 
be very strong when the number of particles is large. The 
shape of the transmitted pulse can differ strongly from 
the initial shape and under certain conditions it gets a 
double peak structure. 

The specific system we consider is a condensate falling 
under gravity and scattering on a Gaussian potential 
barrier. In a recent experiment, it has been shown the 



bouncing of condensates on the mirror formed by a far- 
detuned sheet of light when the strength of the energy 
barrier is much larger than the kinetic energy of the con- 
densate [11]. Here we reproduce the bouncing effect and 
also predict that, when the energy barrier becomes com- 
parable to the kinetic energy, macroscopic quantum tun- 
neling can be observed. We consider the case in which the 
trapping potential in the horizontal directions remains 
active during the fall as well as the case in which the 
trapping potential is turned off at an initial time in all 
directions. We make comparison also with the behavior 
of a one dimensional case. 



II. NUMERICAL PROCEDURE 

The zero-temperature dynamics of a low-density Bose 
condensate can be accurately described by the time- 
dependent Gross-Pitaevskii equation (TDGPE) [12,13] 
for the macroscopic wavefunction (order parameter) 
ip(r, t) of the condensate in an external potential V ex t(r). 
When the wavefunction is normalized to unity then the 
parameter of the nonlinear term of the TDGPE is 7 = 
Airh 2 a s {N — l)/m, where a s the s-wave scattering length 
and N is the number of condensed atoms. 

In our problem, a Bose condensate is initially trapped 
in a harmonic potential and is under the action of gravity. 
In addition there is a Gaussian barrier centered at z = 0. 
Therefore at t < the external potential reads 

V< t {p,z) = -(uy + ^(z - z ) 2 ) + mgz + Ue~^ , 

(1) 

where p = (x 2 + y 2 ) 1 ^ 2 and z are the cylindrical coor- 
dinates, uj p and uj z are the frequencies of the trapping 
harmonic potential, zq is the position of the trap mini- 
mum along the z axis, g is the gravity acceleration, U 
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is the height of the potential barrier and a its width. 
First we consider the case in which an initially isotropic 



trap (Up 



lo z ) is completely removed (to p 



0) 



for t > 0. Then, at positive time, V^, t only contains the 
gravitational term and the barrier. In cylindrical coordi- 
nates, the TDGPE equation becomes 
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which is a nonlinear parabolic partial differential equa- 
tion in 2 + 1 dimensions and we study the time evolution 
starting from the ground state of the system with poten- 
tial (1). 

The numerical integration of this cylindrical TDGPE 
is performed using a modified split operator technique, 
adapted to the integration of a Schrodinger equation [14]. 
We write the Eq. (2) in the form 

d r ~ i 

ifi—i){p,z,t)= yH p (p,z,t) + H z (p,z,t)\ip(p,z,t) , (3) 

where 

H p {p, z, t)=f p + V> t (p, z = 0) + §7|^(p, z, t)\ 2 , (4) 



H z (p,z,t) = T z 



+ V c > xt (p- 



0,z) + -^(p,z,t)\ 2 , (5) 



with T p and T. 



the transverse and axial kinetic opera- 
tors, respectively. The full Hamiltonian is split in two 
sub-Hamiltonians, so that at each time we have to write 
the Laplacian and the external potential with respect to 
one coordinate only, leading to the solution of a tridiag- 
onal system, and to huge savings in computer memory. 
Equation (3) is integrated using the scheme 



iP(p,z,t + 5)= l + A z (t) l-A p {t) 



l + A p (t) l-A z {t) i>{p,z,t) 



(6) 



where S is the integration time step and Aj(t) = 
i8Hj{p,z,t)/(2h) with j = p,z. The splitting is car- 
ried out so that the commutators are exact up to the 
order S 2 included. The nonlinear term j\tp(p, z,t)\ 2 gives 
a problem, because we should really use a tp somehow 
averaged over the time step S, not a evaluated at the 
beginning of the time step. To circumvent this prob- 
lem, we used a predictor corrector method. Each in- 
tegration step t — > t + 5 is performed twice: the first 
time we use ip{p,z,t) in the nonlinear term, obtain- 
ing a predicted ip{p,z,t + 5); we then repeat the in- 
tegration step, starting again from ip(p,z,t), but using 



\ (ipip, z, t) + ^(p, z,t + 5)j in the nonlinear term. In 

this way the solution is accurate to 0(S 2 ). At each time 
step the matrix elements entering the Hamiltonian are 
evaluated by means of finite-difference approximants us- 
ing a typical mesh of 400 x 2400 points in the (p, z) plane. 
As a check of the accuracy of the algorithm we find that 
in the time interval At = IOoj^ 1 the normalization of 
ip is conserved within 5 °/oo and the energy within few 
%. Note that Eq. (3), now with the full potential (1), 
can be used with an imaginary time t = it to find the 
ground-state of the system. 

In our calculations we adopt the harmonic oscillator 
units chosing u H = (c^w z ) 1/ ' 3 = 2ir x 100 Hz. For 23 Na 

atoms, the harmonic length is an = (fi/ '(itiujh)) 1 ^ 2 = 27 
pm and the scattering length is a s = 3 nm. The nonlinear 
term in scaled units is given by 47ra s (iV — 1 )/a#. 



III. RESULTS AND DISCUSSION 

First we study the case in which the initial trap is 
isotropic (uj p = uu z = uor)- In our computations we set 
z = 15an so that the condensate is initially far from 
the Gaussian potential barrier whose effect is negligible. 
Actually, also the gravitational effect is negligible and the 
density profile of the ground-state has a Gaussian shape 
for a small number of atoms and it is an inverted parabola 
for a large number of particle. Then, wc switch off the 
harmonic potential and use the previous wave- function as 
initial condition for Eq. (3). The total energy per particle 
of the condensate is about 180 Hujh with N ranging from 
1 to 10 5 . 
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FIG. 1. Density profile \ip(p = 0,z)| 2 at t = 0, t = 2 and 
t — 4. Barrier parameters: U = 200 and a = 1. iV = 10° 
condensed atoms. Lenght in units an = 27 fim, time in units 




1.6 ms and energy in units Hojh- 



To simulate an impenetrable wall, we have chosen the 
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following values for the parameters of the Gaussian po- 
tential barrier: U = 1000 TujJh and a — an- The center 
of mass motion of the condensate is quasi-periodic due 
to the bouncing of the condensate. Because the conden- 
sate is also expanding isotropically, the oscillation along 
the z axis is damped and this damping increases with 
the number of atoms (larger repulsive interaction). We 
have verified that interference fringes appear when the 
condensate collides on the barrier. 

The phenomenon of interference clearly appears also 
when the energy U of the barrier is comparable to the 
chemical potential /i of the condensate defined as the to- 
tal energy per particle plus the interaction energy per 
particle. The chemical potential, in contrast to the en- 
ergy, is time-dependent [15]. Nevertheless, in our calcu- 
lations fi is practically constant, apart a slight increase 
(no more than 4%) during the collision time due to the 
enhanced density of the condensate. In Fig. 1, we plot 
the density profile of the condensate along the symmetry 
axis at three different instants. In this case the energy 
barrier is U = 200 huj H and the initial chemical potential 
fj, = 184.23 Uuj h . Besides the interference between the in- 
cident wave- function and the reflected wave- function, one 
sees a transmitted wave, due to the quantum tunneling 
of the condensate. 

For a non-interacting gas the scmiclassical ap- 
proximation of the linear Schrodinger equation pre- 
dicts that the tunneling probability is proportional to 

e 2 ^ dz V 2m ( v ext( z ) ; w here a and b are the clas- 

sical turning points and E is the energy per particle of 
the condensate that, in the non-interacting case, coin- 
cides with the chemical potential fi. In the interacting 
case, on the basis of a WKB study of the ID GPE in the 
weak-coupling regime [8] , it has been shown that the en- 
ergy per particle E should be substituted by the chemical 
potential \i of the Bose condensate. 

In general, the semiclassical quantization of the non- 
linear Schrodinger equation appears a difficult task but 
numerical results can be easly obtained with our algo- 
rithm. As shown in Fig. 2, the effect of nonlinear term 
is such that the tunneling fraction Pt increases with the 
number N of condensed atoms, i.e. by increasing the 
chemical potential. For instance, for N — 10 5 the tun- 
neling probability is about 4 time larger than the non- 
interacting case. Note that, due to the high density in the 
impact region, the self-interaction becomes quite strong 
for a large number of particles and this clearly affects the 
tunneling fraction, as shown in Fig. 2. 

While the transmitted wave- function of the condensate 
falls towards z = — oo, the reflected wave-function goes 
back close to the initial position zq and then falls down 
again under the force of gravity. It means that, after a 
time roughly equal to the period 2y / 2za/g, the conden- 
sate scatters again on the Gaussian barrier and now the 
tunneling fraction is smaller. 



0.1 



0.08 - 

N=10 5 

N=10 4 

0.06 - N=10 3 

N=1 

0.04 - 



0.02 - 








12 3 4 

t 

FIG. 2. Tunneling fraction Pt as a function of time t. 
Barrier parameters: U = 200 and a = 1. N — 1 means 
non-interacting case. Initial position: zq — 15. Units as in 
Fig 1. 

In Fig. 3 we plot the tunneling fraction Pt as a func- 
tion of time, showing that the macroscopic quantum tun- 
neling of a Bose condensate falling under gravity is a 
quasi-periodic phenomenon. This shows that experimen- 
tal set up similar to what we are considering should be 
able to generate coherent atomic pulses. 
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FIG. 3. Tunneling fraction Pt as a function of time t. Com- 
parison between radially-free (lo p — 0) and under radial con- 
finament (co p — 0.5) falling condensates. N = 10 5 condensed 
atoms and initial position zo = 15. Barrier width: a = 1. 
Units as in Fig 1. 

In the previous computation the condensate expands 
for t > in all directions but this leads to a decrease in 
density at the second bounce with a reduced tunneling 
fraction. It is interesting to see what happens if the con- 
densate remains trapped in the transverse directions. We 
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obtain this assuming that V^. t does not confine (u) z = 0) 
in the z direction but it remains confining (uj p = auin) 
in the horizontal directions. 
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FIG. 4. Tunneling fraction Pt as a function of the ini- 
tial chemical potential fx for different frequencies lo p of ra- 
dial confinanicnt. The points in each line correspond to 
a sequence of numbers N of atoms. From left to right: 
N = 1,10 2 ,10 3 ,10 4 ,5 x 10 4 ,10 5 . Barrier parameters: 
U = 200 and a = 1. Units as in Fig 1. 

From Fig. 4 we see that, depending on the degree of 
this radial confinement, the interatomic interaction (and 
the chemical potential jj,) either enhances or depresses 
the tunneling probability Pt- Note that, in the present 
case, the system is not quasi- ID, due to the isotropic ini- 
tial condition. In the case of strong radial confinanicnt, 
the falling condensate shrinks in the radial direction dur- 
ing the motion. Such compression modifies the effect of 
the nonlinear term in the TDGPE. As shown in Fig. 4, 
the tunneling probability Pt initially decreases and then 
increases with the chemical potential ll. 

In the previous calculations we have considered Bose 
condensates with isotropic initial conditions (uj z — u> p ). 
Often experimentally the confining potential of the trap 
is anisotropic. Here we investigate the case of cigar- 
shaped initial conditions (u> z = u) p /10) which model a 
quasi-lD system. In this case we find that the tunneling 
fraction increases by increasing the number N of parti- 
cles both in the case of absence of confinement (u> p = 0) 
and in the case of strong radial confinement (Fig. 5). As 
shown in Fig. 5, due to the anisotropy of the trap, the 
axial elongation of the condensate grows with the number 
N of particles and for N = 10 4 one tip of the condensate 
touches the Gaussian barrier. 
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FIG. 5. Density profile \ip(p = 0, z)\ 2 at t = for various 
values of the number N of condensed atoms. Barrier param- 
eters: U — 200 and a = 1. Radially-confined (u) p — 1) falling 
condensate with cigar-shaped (lu p /lu z = 10) initial condition. 
Pt is the tunneling fraction. Units as in Fig 1. 

The exponential law suggested by Zapata et al. [8] for 
a ID system seems able to capture the behavior of Pt 
as a function of the chemical potential li: on the basis of 
this law one finds a relative increment of the tunneling 
fraction of 87% as N changes from 1 to 10 4 , whereas 
this relative increment is about 81% on the basis of the 
TDGPE. 

In Fig. 6 we plot the density profile of the condensate 
at three time step. One sees that the transmitted pulse 
has a peculiar feature: the density profile clearly shows 
two peaks (see Fig. 6). The amplitude of the second 
slower peak increases with the number of particles. 
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FIG. 6. Density profile \ip{p = 0,z)\ 2 at t = 0, t = 2 and 
t — 4. Barrier parameters: U = 200 and a = 1. iV = 10 3 con- 
densed atoms. Radially-confined (cjp = 1) falling condensate 
with cigar-shaped {u) p /oj z = 10) initial condition. Units as in 
Fig 1. 

Actually, for a large number of particles a third small 
peak can be detected. In the case of an isotropic cloud, 
the second peak shows up only with N > 10 5 particles 
(see Fig. 1). These two features of the tunneling of a 
cigar-shaped cloud, a tunneling rate enhanced by the in- 
teraction and the presence of a second peak in the trans- 
mitted pulse, are also present in the case of a purely one 
dimensional system. 

For the sake of completeness, we have also considered 
the single-vortex state ^>(r, t) — p<p(p, z)e %e as initial con- 
dition but we have not found any significant difference in 
the dynamics of the condensate. 

We have presented the results for the Na atom. As long 
as the scattering length is positive our results apply also 
to the other alkali-metal atoms by suitably rescaling the 
spatial and temporal scales and the number of particles. 



IV. CONCLUSIONS 

We have studied the dynamics of Bose condensates 
falling under gravity and scattering on a Gaussian barrier 
that models a mirror of light. Apart the pure bouncing 
with interference, that can be seen for very large val- 
ues of the energy barrier, we have investigated quantum 
tunneling. Our results show that macroscopic quantum 
tunneling of a Bose condensate should be observable ex- 
perimentally and that the tunneling rate is very sensitive 
to the interatomic interaction and to geometrical features 
of the experiment. If initially the cloud is spherical under 
free fall, the tunneling fraction grows with the number 
of particles, which increases the interatomic interaction 
and the chemical potential of the condensate. But, the 
fraction is reduced by the interatomic interaction if the 
falling condensate is under a strong enough transverse 
confinement. In a ID case, and in the case of a conden- 
sate with a cigar-shaped initial condition, the interatomic 
interaction always enhances the tunneling rate. There- 
fore the tunneling rate becomes larger and larger as the 
number of particles increases. In addition, in the ID or 
quasi ID case the density profile of the transmitted wave 
shows two peaks. We have also shown that in our system 
macroscopic quantum tunneling is a quasi-periodic phe- 
nomenon and it can be used to generate Bose condensed 
atomic pulses. One could generate truely periodic pulses 
if the bouncing part of the condensate is reconfined for a 
short period of time when it returns in the initial position 
and additional atoms are condensed in order to regener- 
ate the initial state. Often the dynamics and static of a 
condensate is modelled by a ID system. One important 
result of our computation is that geometrical aspects are 



important in tunneling phenomena and that ID results 
not always reflect the behavior of 3D systems. 
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